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T-H Abstract 

o 

^vj We characterize the small-time asymptotic behavior of the exit probability of a Levy process out of a two- 

sided interval and of the law of its overshoot, conditionally on the terminal value of the process. The asymptotic 
expansions are given in the form of a first order term and a precise computable error bound. As an important 
application of these formulas, we develop a novel adaptive discretization scheme for the Monte Carlo computation 
of functionals of killed Levy processes with controlled bias. The considered functionals appear in several domains 
of mathematical finance (e.g. structural credit risk models, pricing of barrier options, and contingent convertible 
bonds) as well as in natural sciences. The proposed algorithm works by adding discretization points sampled from 

n^ the Levy bridge density to the skeleton of the process until the overall error for a given trajectory becomes smaller 

^ , than the maximum tolerance given by the user. 

r^^ Keywords and phrases: Small-time asymptotics. Levy bridge, killed Levy process, exit probability, bridge Monte 

jrt Carlo methods, adaptive discretization, barrier options. 
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^ 1 Introduction 

'^ Small-time asymptotics for the distributions of Levy processes and related Markov processes have a long history going 
^i back to the seminal work of Leandre [30l , who obtained the leading order term of the transition density of a Markov 
process solving a stochastic differential equation with jumps. In the case of a Levy process, the main result of [30j 
^2J reads 
^ lim^Mx) = six), (x^O), (1.1) 

tj where ftix) :— j^F{Xt < x) is the marginal density of the Levy process X and s is the Levy density of X, whose 

. ^H existence and smoothness need to be assumed. Leandre's approach was to consider separately the small jumps (say, 

^v those with sizes smaller than an e > 0) and the large jumps of the underlying Levy process, and to condition on the 

J-j number of large jumps by time t. A similar approach has been applied during the last decade to obtain high-order 

asymptotic expansions for the transition distributions and densities of Levy processes and other Markov processes with 

jumps (see [38], [19], [20], and [21]). These small-time asymptotic results have found a wide scope of applications 

ranging from estimation methods based on high-frequency sampling observations of the process (see, e.g., [17], [TT] . 

[37], and references therein) to asymptotic results for option prices and Black-Scholes volatilities in short-time (c.f. 

[S], [IS], [IS]). 

In the present paper, we adopt Leandre's approach to study the asymptotic behavior of the generalized moments 
of the Levy process stopped at the time it exits a two-sided interval (a, b), conditionally on the terminal value of the 
process. Specifically, for a Levy process {Xt)t>o with Levy density s that is smooth outside any neighborhood of the 
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origin and for a bounded Lipschitz function (/s, we prove that 

E{^{X,)lr<t\Xt = y) = l [ ^(«)^^^^lfc^dz; + o(t), {t ^ 0,y e {a,b)\{0}), (1.2) 

2 7(„b)c s{y) 



where r := inf{u > : X„ ^ (a,^)} with — oo < a < < 6 < oo. In the case tp= I, (1.2) can be written as fohows 



¥i3ue[0,t]: X^(^{a,b)\Xt = y) = l f '-^"^^^ dv + o{t) , (t^O), (1.3) 

for y e (a, 6)\{0}. As in the case of the small-time asymptotics for the marginal distributions of the process, the main 



intuition can be drawn from considering the pure-jump case with finite jump activity. Intuitively, formulas (1.2)-(1.3) 
tell us that if, within a small time period, a Levy process goes out of the interval (a, b) and then comes back to the 
point y € (a, &), this essentially happens with two large jumps: the first jump takes the process out of (a, b), while the 
second jump brings it back to y. 



Our study of the short-time behavior of (1.2) and (1.3 1 is motivated by applications in the Monte Carlo evaluation 
of functionals of the form 

E[F{XT)lr>Th T = inf{t > : X, ^ (a, b)}. (1.4) 

In financial mathematics, such functionals arise in structural credit risk models based on Levy processes |16j and in 
the pricing of barrier options (cf. [27], [7]), which is one of the most popular classes of exotic options. Very recently, 
a renewed interest to these problems has emerged in relation to the so-called contingent convertible bonds, where the 
conversion is triggered by a passage across a level and which exhibit a high sensitivity to jump risk |13j . In natural 
sciences. Levy processes (under the name of Levy flights) are used as models for certain diffusion-like phenomena in 
physics and chemistry (so-called anomalous or super-diffusion) |32l 1411 13] as well as to describe movement patterns 
of foraging animals [44l [5] , and there is considerable interest towards the study of Levy flights in bounded domains 
and related first passage problems giving rise to functionals of type ( |1.4[ ) [101 [HI HI]- In all these settings, closed- form 
expressions are rarely available and Monte Carlo is often the method of choice. 



The simplest procedure to evaluate the functional (1.4) by Monte Carlo consists in simulating the process (Xt) 



t>o 

at evenly spaced times ij! := fc/i„, with /i„ :— T/n and k = 0, . . . , n, over the interval [0, T], and approximating the 
exit time t by 

f„:=inf{t^X,. ^(a,5)}. 

This simple method introduces two types of errors: the statistical error and the discretization error. The latter is 
known to be quite significant (cf. [2] and Example 2 in Section [5] below); [31 reports errors of up to 10% in the context 
of barrier options for a time discretization of one point per day. 

In the context of continuous diffusions, short-time asymptotics have been successfully employed to alleviate the 
bias due to the discretization error. One of the earliest procedures of this type, due to Baldi [5], is based on an 
approximation of the probability, p(x, y, t), that the process X has gone out of a domain (a, b) during the small time 
interval [s, s -I- i] conditioning on Xg = x and X^^t = y] i-©-, 

p(x, y,t):^Vi3ue[s,s + t]: Xu(^ia,b)\Xs= x, X^+t = y) ■ (1.5) 

Given such an approximation p{x,y,t) of the functional p(x,y,t), the procedure simulates iteratively Xfn at each 
step A: = 0, . . . , n — 1, and if Xt^^ G (a, b), it proceeds to kill the process with probability p{Xt^ , Xj" , /i„) and choose 
f^j^i = (fc -f l)hn as an approximation of the exit time r. A similar idea was used in [33^ to price barrier options with 
payoff (f{Sr,T) by Monte Carlo. 

In the context of Levy processes, an attempt to apply a similar methodology has been made in [?S1[3S]. The authors 
remarked that the discretization bias can be reduced by using the identity 



and replacing the exact exit probability p{x,y,t) with a suitable small-time approximation p{x,y,t). However, these 
papers propose no general formula for p{x,y,t) and, as shown in [3], the Monte Carlo method proposed in [451 136) 
could lead to a large discretization bias. On the other hand, in the specific case of the parametric variance gamma 
model, there exist discretization algorithms (cf. [1]) allowing to simulate the running minimum and maximum with 
error bounds. Let us also remark the recent work of [28j where a method for the joint simulation of the running 
maximum and the position of a Levy process is introduced based on the Wiener-Hopf decomposition of the process. 



Our short-time asymptotic result (1.3) provides an approximation of the exit probability (1.51 via the formula 



p{x,y,t) 



s(v)s{y — X ~ v) t 

r aw = - 

(a-xA-xY s{y-x) 2 



{a,bY 



s{u — x)s{y — u) 
s{y-x) 



du. 



(x^y), 



(1.7) 



which is valid under mild regularity conditions on the Levy process X (see Section [2] for details). The first order 
approximation (1.7), together with an appropriate error bound for it, enable us to develop a general adaptive Monte 



Carlo method for evaluating the functional (1.4) with a given precision. Given a target error level 7, the idea is to 
generate a "random skeleton" {{Tk,XTk)}k=i of the process X such that the error in each subinterval [Tk,Tk+i], i.e. 



■ Pi^n 1 ^Tfc+i , Tfc+i -Tk) - p{Xt^ , Xt^^^ , Tfc+i - Tk) 



satisfies |e| < -^^ — ^7. The functional (1.4) is then approximated as follows: 



N~l 



nF{XT)lr>T]^^\F{XT) n {l-p{XT,,XT,^,,n+l-n)] 



(1.8) 



(1.9) 



fc=0 



and it is shown that the total bias of this computation will be less then 7. As a result of this adaptiveness, the algorithm 
generates more frequent points when the process X is close to the boundary, and takes large time steps (thus saving 



computational time) when the process is far from the boundary. Let us re mark that, unlike the formula (1.6), where 
the sampling times {t^} are deterministic and fixed, the decomposition ( 1.9 ) for random skeletons X :— {{T^, XTt,)}^=o 



requires precise (and also novel to the best of our knowledge) conditions under which this formula holds (see Section 
|4]for the details). 

The proposed adaptive algorithm works as follows. First, the endpoint X^ is generated and added to the skeleton. 
Next, if the error (1.8) is too large for a given subinterval [T^jT^+i], the procedure splits the interval into two and 
generates the midpoint Xf^ with Tk :— {Tk + Tk+i)/2 from the bridge distribution. This is repeated iteratively until 
the desired error bound is satisfied for every subinterval [Tk,Tk+i] of the sampling times = Tq < ■ ■ ■ < T/v = T. 
Such retrospective sampling (starting from the endpoint) has a number of advantages over the classical uniform 
discretization, especially in the context of rare event simulation, where it enables one to easily implement variance 
reduction by importance sampling. Indeed, the process can be directed to the region of interest by modifying the 
distribution of the terminal value, while keeping unchanged the rest of the algorithm. On the other hand, this method 
requires fast simulation from the bridge distribution oi Xt/2 conditioned to Xt = y. To this end, as another contribution 
of particular interest on its own, we also propose a new method to simulate from this Levy bridge distribution based 
on the classical rejection method. 

As previously explained, in order to implement the above adaptive algorithm, precise computable bounds for the 



approximation errors in (1.2)-(1.3) are also needed. We obtain such bounds by developing explicit inequalities for the 



tail probabilities and transition densities of a Levy process whose Levy density has a small compact support. This type 
of concentration inequalities in turn allows us to estimate the different components of the error, which, as explained 
above, originate from conditioning the desired functional on the number of big jumps by time t (see Section p^ for the 
details). The resulting error bounds are given in terms of the Lipschitz and Loo norms of ip as well as several computable 
quantities related to the Levy density s such as sup|2,|>g s(x), supu|>£ |s'(x)|, /|^|>g s{x)dx, and /|^|<j x^s{x)dx. 

Let us also remark that an adaptive simulation method similar to the one introduced in the present paper was 
proposed in [15] to compute a functional of the form Kip{Xt-, t) for a homogeneous diffusion process X without jumps. 
Adaptive numerical methods for finding weak approximation of diffusions without jumps and with finite intensity jumps 
(but with the adaptiveness only concerning the diffusion part) have also been proposed in [42] and |34j, respectively. 
As in our paper, the idea therein is to sample from inside of a subinterval [i}J,i}J_]_i] whenever the approximation error 
in that subinterval has not reached a desired low level, specified by the user. 



The paper is organized as follows. In Section[2l we obtain the leading term of the functional E {Lp{XT-)lT-<t\ Xt = y) 
when i -^ 0. The explicit estimate of the approximation error is given in Section |3] The development of the adaptive 
discretization schemes for the Monte Carlo computation of the functional ¥\F{Xt)'^t>t\ as well as the algorithm to 
simulate random observations from the Levy bridge distribution are given in Section |4J Our methods are illustrated 
numerically in Section [5] for Cauchy process. Finally, the proofs of the technical results are deferred to the appendix. 

2 Small-time asymptotics for Levy bridges 

Let X be a real- valued Levy process on a probability space {Vt^F.V) with Levy triplet ((T^,j/, /i) with respect to 
truncation function h[x) = l|a;|<i. Throughout, {Tt)t>o denotes the natural filtration generated by the process X and 
augmented by the null sets of J-" so that it satisfies the usual conditions (see, e.g.. Chapter L4 in f35^). The following 
standing assumptions are imposed throughout the paper: 

• The Levy measure v admits a continuously differentiable density s : ]R\{0} -^ (0, oo), with respect to the Lebesgue 
measure (hereafter denoted by £), which satisfies 

sup s{x) < oo, sup |s'(x)| < oo, (2.1) 

for any e > 0. 

• The distribution of Xf admits a density ft for all i > 0. Since v is already assumed to admit a density, for this 
assumption to hold, it suffices to additionally require that j/(K) = oo or cr > (see Theorem 27.7 in Sato |1D]). 

• The density of Xt satisfies ft{x) > for all a; £ M and t > (see Theorem 24.10 in Sato [JOj for mild sufficient 
conditions for this property to hold). 

As it is usually done with Levy processes, we shall decompose X into a compound Poisson process and a process 
with bounded jumps. More specifically, for any e e (0,1), we select a function c^ € C°°(M), which is decreasing on 
(— oo, 0) and increasing on (0, oo) and such that l|a;|>e < Ce(a;) < l\x\>e/2- Next, we define the truncated Levy densities 

Se{x) := Ce{x)s{x) and Se{x) := Ce{x)s{x), 

with Ce{x) := 1 — c^{x). Let Z^ be a compound Poisson process with Levy measure Ss{x)dx and X"^ be a Levy process, 
independent from Z^ , with characteristic triplet {<T'^,s^{x)dx, ^e), where 



11^ := fi — / xcs{x)s{x)dx. (2.2) 

"'|2:|<l 

It is clear that X^ + Z'^ has the same law as X and that the intensity and probability density of the jumps of Z^ are 
Ae :~ J Se{x)dx and Ss{x)/Xe, respectively. Throughout the paper, we let {Nt)t>a be the jump counting process of 

Z^ and {Y^)k>i be the jump sizes of Z^. Thus, Zf — X]/e=i ^k ■ Note that the distribution of Xf is also absolutely 
continuous since cr > or J s^{x)dx = oo, for any e > 0. For future reference, let us remark that 



E {Xf) = t [ ^j + / xs^{x)dx j = i^£, Var {Xf) =ticj'^+ / x'^s^(x)dx\ =: ta 



(2.3) 



\x\>l 

since e S (0, 1) (see, e.g.. Example 25.12 in [10] for the mean and variance formulas of a Levy process). 

The following Lemma will be needed in what follows (c.f. Propositions 1.4 and III. 2 in [5IP). See also Sections 
3.1|3.2 below for explicit expressions for the constants Cp{rj,£) and Cp(77,e). 



Lemma 2.1. Let ft he the transition density of the small-jump component process {Xl)t>o- Then, for any fixed positive 
real rj and positive integer p, there exist an eo{ri,p) > and positive constants tQ(ri,e), Cp{ri,e), and Cp{rj^e) < oo for 
any e < Eq such that 



(i) P sup \XI\ > ry < Cp{7j,e)tP, (ii) sup f!ix) < Cp{7j,e)tP, (2.4) 

V0<s<t / |a:|>i) 



for all < t < to and < e < Eq. 



The following result provides the key tool for establishing the small-time asymptotics of the moments of the Levy 
bridge "stopped" at the exit time from an interval (a, b). Its proof is presented in Appendix [A| 

Theorem 2.2. For fixed constants a e [— cx),0) and b G (0,cx)], define 

r :=inf{u> : X„ ^ (a,b)} . 

Let c^ : M — >■ M be bounded and Lipschitz on M. and let Sq £ (O, — ^) • Then, for any y £ (a + (5o, 6 — Sq) and < S < Sq, 



ry+S / £2 r 

E {ip{Xr)l{r<t,Xte{y-5,y+S)}) = [^ 

Jy-S \ ^ J(a,l 



ip{v)s{v)s{u — v)dv + TZt{u)t^ du, 



where the remainder term TZt{u) is such that 



lim ess sup \TZt{u)\ — 0. 



t-i-O 



ue(a+5o,b—So) 



(2.5) 



(2.6) 



Remark 2.3. By the definition of conditional expectation, 

E {ip{Xr)l{r<t,Xte{y-5,y+S}}) ^ E (E {ip{Xr)l{r<t}\ ^t) i{Xte{y-S,y+S)}) 



y+s 



y-5 



E {(p{Xr)l{r<t}\ Xt = u) ftiu)du, 



(2.7) 



where ft{u) is the density of Xt and, as usual, $(w) := E ((p{XT-)lt^^n\ Xt — u) is such that ^{Xf) is a version of 
E (^(fi{XT-)'i-{r<t}\ ^t) ■ Comparing (2.1) and (2.5), it then follows that, for C-a.e. yE {a,b). 



E{^{Xr)l{r<t}\Xt^y) 



Tl{a,br Vi'")-'^iv)s{y - v)dv nt{y)t^ 



ft{y) ft{y) 

If, in addition, the transition density ft satisfies the asymptotic formula (1.1 ^then, for C-a.e. y £ (a, 6)\{0}, 



[, ,,^ u:>(v)s(v)s(y — v)dv 
E{^{X^)^^^t}\Xt = y)^t ^^^^'^''^' '^'^^^^'" +o{t). 



(2.8) 



(2.9) 



Formulas \2.b\ and (2.8 I can be interpreted as large deviation results for the trajectories of Levy processes in small 



time. When ip{x) = 1, (2.9) gives the following small-time approximation for the exit probability of the Levy bridge. 



'{T<t\Xt = y)=t 



2s(2/) 



o{t). 



(2.10) 



We conclude this section with a simpler result for the case when Xt is outside the interval. Its proof is outlined in 
Appendix |Xj 

Proposition 2.4. Let y) : M — >■ K be bounded and Lipschitz on M, and let 5q > 0. Then, under the same notation and 



conditions as in Theorem 2.2 for any y € {a — So,b -\- SqY and S < 6o, 



ry+s 
E (ip{Xr)l{Xteiy-s.y+S)}) = / {tip{u)s{u) + nt{u)t) du, 

Jy-s 



where the remainder term TZt{u) is such that 



lim ess sup |''^t(u)| = 0. 



t->o 



uefa—So^b+So)" 



(2.11) 



(2.12) 



Remark 2.5. Analogously to Remark 2.3, (2.11) enables us to establish the following natural asymptotic formula: 



E{^{Xr)\Xt ^y)^t ^^M^ + o(l) - ^(y) + o(l), {t ^ 0), 

.ft{y) 

for C-a.e. y E [a,bY. The second equality above holds whenever ft{y) satisfies (UTTp. 

^As stated in the introduction, | |l.l[ | liolds for a large class of Markov processes with jumps as proved by 1301 . For Levy processes, 1381 
provided a more elementary proof using the same conditions and similar approach as in I30| . Higher order short-time expansions for the 
transition densities were obtained in |19| . 



3 On a precise bound for the remainder term 

In the previous section, we developed the necessary resuhs for finding estimates of the functional 

f{0,y,t) ■.= E[ip{Xr)lr<t\Xt ^ y] (3.1) 



in short-time. Indeed, as explained in Remark 2.3 Theorem 2.2 yields the following natural estimate for f(0,y,t): 



ft[y) 

The estimate ( |3.2[ ) will be used below to develop adaptive discretization schemes for the Monte Carlo computation of 
functionals of the killed Levy process (see Section l4|. To this end, we first need to find an explicit estimate for the 



remainder TZtiv) appearing in (2.5 1. Such an estimate can be expressed in terms of bounds for the tail probability and 



transition densities of the small-jump component {X[)t>o- Hence, we start by providing explicit expressions for the 



upper bounds appearing in (2.4) and then proceed to give a precise error bound for \f{0,y,t) — /(O, j/,t)|. 



3.1 Bounding the tail probability of the supremum 

The following exponential inequality for Levy processes with bounded jumps will be important to estimate the supre- 
mum of the small-jump component (X^) defined in Sectional Its proof, which is provided in Appendix [B] for com- 
pleteness, is a variation of the bound obtained in [38j (which in turn is based on Lemma 26.4 in j40| ) . 

Lemma 3.1. Let {Mt) be a martingale Levy process with lAM^j < e and (Af, Af)t = a^t. Then, 

sup(Af, +Ais)>77) <i?C,(r;,e;//) (£ = 0,1), (3.3) 

\s<t / 

with the following constants Ci{ri^e]^) and corresponding conditions: 

(1) Co{ri,e;fi) = e^^ e"^^'^ for all r/ > and < t < r]/[^ V 0) (with the convention here and below that the 
fraction is +cx) if the denominator is zero); 

(2) Ci{T],e;n) = e^" ' (^) for all i] > and <t < ri/{^\JQ) if either (i) fi < or (ii) ^ > and ij < al/e; 



In order to apply Lemma 3.1 for (Xf)t>o, we recall that < e < 1 so that EXf = /i^t. Then, the martingale part 



Mf := X^ - n^t of X^ is such that 

(A'r,M'^)j = {a^+ I c,{x)x^i^{dx) ] t = a^t. 
Thus, fixing 



2(Me V 0) 



it follows that, for all < i < toi 



'supXf >v)<¥ (supAf« + \l^e\t >v]<P (supA^f > i) < **^ (1'^) 
with C{ri, e) is defined by 



(3.5) 



C(ry,.):=(^) . (3.6) 



eat 



Similarly, we have P(sup,.<t \Xl\ > ?;) < 2iAC (?7/2,£). 



3.2 Bounding the transition density of the smah-jump component 



To obtain explicit expressions for the constants appearing in the bounds for the density f^ in Lemma |2.1[ we shaU 
assume that the process X is such that X^ has a unimodal distribution for ah i > and £ > 0. By Yamazato's theorem 
(see Theorem 53.1 in [1^), a sufficient condition for this is that the process X is self-decomposable, which is the case if 
and only if the Levy density s is of the form s{x) = -j^ for a function k which is increasing on {— c», 0) and decreasing 
on (0,oo) (see Corollary 15.11 in [40]). In particular, most of the parametric models used in the literature (such as 
stable, tempered stable, variance gamma, and normal inverse Gaussian processes) are self-decomposable and so these 
processes as well as their truncated versions have unimodal densities at all times. 

Let mf be the mode of Xf . If mf G [—77, 77] and r; > 77, then the density can be estimated by 



sup f^{x) < 

\x\>r, V-IJ 



P[\x!\>vi 



(3.7) 



simply because the density is decreasing in (77,00) and increasing in (—00,-77). The relation (3.7) in turn leads to a 
bound of the form (2.4-ii) by applying the tail bound (2.4i). It remains to find conditions for m^ € [—77,77]. Since 
obviously X^ has finite second moment, the following bound due to Johnson and Rogers [2& can be applied 



\ml-EX',\'<3Y^T{X^). 



(3.8) 



Thus, recalling the mean and variance formulas given in (2.31, mf G [^VtV] whenever < i < ii, where ii is such that 

ti\ti,\ + V3tY^ la^+ [ \x\^i^{dx)] =77. (3.9) 

\ JN<e J 

By taking 77 = 77/2, we will have 

4 r m fin(n/A r) „ 

(3.10) 



sup ft{x) < -1 

\x\>r, V 



\xn > I 



^ 8C(77/4,£) ^^^ 



for any < t < ti A to with to defined as in (3.4) 



3.3 Precise bound for the remainder 



We are now ready to give an explicit bound for the reminder term TZt{y) appearing in ( 2.5 ), which in turn will produce 



an error bound for |/(0, y, t) — /(O, y, t)\. Throughout, we shall use the following notation: 

(i) Oe := supj. Se{x) and a'^ := sup^ |sg(a;)|, where, as before, Se{x) := Ce{x)s{x) is the Levy density s, truncated in 
a neighborhood of the origin; 

(ii) Ag := J s{x)ce{x)dx, 11^ := ^ — S\x\<i xce{x)s{x)dx, and a^ := cr^ + J c^{x)x'^ s{x)dx; 



(iii) C{r],e) is defined as in (3.6), to{s,r]) is defined as in (3.4), and ti{e,T]) is defined as in (3.9) 



The following result, whose proof is given in Appendix [B[ gives an estimate for TZt{y) in terms of the previously 
defined notation and the Lao- and Lipschitz norms of (p denoted hereafter by 



m\ 



:=esssup|(^(x)|, \\Lp\\up ■= ini {K > : \ip{x) - ip{y)\ < K\x - y\, yx,yeR}. 



Theorem 3.2. Using the notation of Theorem 2.2, assume that the process X is sueh that X^ has a unimodal 
distribution for all t > and e > 0. Let c := b A \a\ and Ay :~ (b — y) A {y ~ a) > 0. Then, 

\'}ltiy)\<^en{0,y,t), 



for alio <t < to{e, {Ay/2) A c) A ti(e, Ay/2), where 



eniO,y.t) := e~^^*M^C{Ay/4,s)t^ 



A„ 



+ 2aet + a^Xer \ + 2e~-^^'||(/7||ooaeC(c/2, £)t^+K {1 + iAJ 



V" + ||</j||ooaeA7i (1 - e-^=*[l + Aei + (Aei)V2]) 



(3.11) 



A„ 



Two immediate conclusions can be drawn. First, note that, by taking e < -^^ A |, we obtain a bound for the 



remainder satisfying condition (2.6). Second, as seen in Remark 2.3 the previous bound imphes the following error 
bound 

1/(0, y, t) - /(O, y,t)\< ^^||yy^ =: e;(0, y, t), 



with / and / defined as in (3.1)-(3.2) 



Remark 3.3. The approximation for the conditional exit probability p{0,y,t) :~ V[T<t\Xt=y] is obtained by 
substituting (p = 1 into (2.8).' 

My) ■ 



p{o,y,t) 

Making this substitution in the previous bound, it follows that \p{0,y,t) — p(0,y,t)\ < ep(0,y,t) with ep(0, y, i) given by 



ep(0,y,i) := 



1 



My) 



£"£+3 



e-^'*C{Ay/A, e)tT^ I — + 2aet + a^X^t^ \ + 2e"^'*aeC(c/2, e)t^+^- {1 + tX^} + -^t 
+ a,X-' (1 - e-^^*[l + A,i + (A,t)V2]) + e-^'H^ [2al + A,<] {a,t^'^ + ^i) 



valid for all t < tQ^s, {Ay/2) A c) A ti(e, Aj,/2). The one-sided case (a — — ooj can similarly be obtained. 



4 Adaptive simulation of killed Levy processes 

Our goal in this section is to design a type of adaptive Monte Carlo estimators for functionals of the form 

E[F{XT)lr>T], 



(4.1) 



where i^ is a Borel measurable function and r :— mi{t > ■ Xt ^ D} with D :— {a,b), for some a G [— oo,0) and 
6 £ (0, oo]. From now on, to simplify notation and with no loss of generality, we shall take T = 1. 

For < s < t, X e M, and y e M, wc denote by PfJl ^ „) [•] the bridge law of the Levy process X on the time interval 
[s,t] with starting value x and terminal value y; that is, this is a version of the regular conditional distribution of 
{x + Xu-s}u£[s,t] given Xt_s = y — x. Since Xt has a strictly positive density on M for every i > 0, the bridge law is 
uniquely defined for £-almost every y G M (recall that C stands for the Lebesgue measure), which is sufficient for our 
purposes. We also let p{x, y, t) denote the exit probability from the domain D before time t for the Levy bridge: 



p{x, y, t) := Ffolccv) [t < t] ^V [3u e [0,t] : x + X„ ^ {a, b)\Xt^y]. 



Our approach is based on the following decomposition: 



E[i^(Xi)l,>i]=E 



F{X,) H {1-p{Xt^,Xt,^„T,+i-T,)) 



(4.2) 



(4.3) 



where = Tq < • • • < Tjv — 1 are suitable sampling times. Formula (4.3) directly follows from the Markov property 
when the sampling points are deterministic. In that case, the set of points X := {{Ti, XTi)}ilo i^ called a deterministic 
skeleton. In our setting both the number of points N and the sampling times = Tq < Ti < • • • < T/y = 1 are random 



and we need to formalize under what conditions on X (4.3 ) still holds. The following result will suffice for our purposes. 



Lemma 4.1. Let N be a random variable with support A/" C N, such that N > 0, and let = Tq < ■ ■ ■ < T/v ^ 1 be 
random points such that 

(1) Each Ti takes values in a countable set K. d [0, 1]; 

(2) For each n £ J\f and {sq, . . . ,Sn) G K-"^^ with — sq < ■ ■ ■ < Sn — ^, the event {N — n, (Tq, . . . ,T„) — 
(sq, . . . , s„)} is o^iXs- : z = 0, . . . , n) -measurable. 



Then, {4-3) is satisfied for any measurable function F with E.[\F{Xi)\] < cx) and, furthermore, for every t G (0,1), 
n€j\f, andAe 6(M), 

¥[Xt£A\N = n,To,...,TN,XTo,--.,XTj=¥^^'}^T,,+„XT,,,XT,,+A^teA], where i* = maxjz : T, < i}- (4.4) 

Proof. Throughout, we let p{x,y,t) := 1 — p{x,y,t), IC^ := {(so,---,Sn) G JC"'^^ : = sq ^ • • ■ < Sn = 1}, 
Uq := [so, si], and Ui := (si, s^+i], with i = 1, . . . , n — 1. We also use the notation 



lu := l{x„e(a,6):Vue;7}, for a domain U C 1R+, and It/) = 1. 
Then, by Markov property 

E[F{X,)lr>l] ^ Y. E E[f(Xi)I[o,i]l{^.„,(T„,...,T„) = (.o,.....„)}] 

"e^(so,...,s„)e/C" 

rn-l 



(4.5) 



E E E 

"e^(si,...,s„)eK" 



-F'(^l)l{Ar=«,(To,...,T„) = (so,...,s„)}E 



n^t/.l^s, :j=0,...,i 



1=0 
n-1 



E E E 

"e^(si,...,s„)e/C" 

Af-l 

E F(Xi ) Y[ p (Xt. , Xt.+, , r,+i - T, 



-F'(-'^l)l{Ar=n,(To,...,T„) = (so,...,s„)} J_ J_ P (-^T, , -'^Ti+i , ^i+l - ^i. 



1=0 



i=0 



which proves (4.3). Similarly, 



F[Xt e ^|A^ = n, To, . . . , Tjv, Xt,,, . . . , Xt„] 

^ P[Xt e A|iV = n, (To, . . . , r„) = (so, ■ . . , s„), Xt„ . . . , ^t„]1(To,....t„)=(.o....,.„) 
(sQ,...,s„)e^" 



(so,.--,sr.)eK:" 



D 



From (4.3), it is now evident that, for the computation of (4.1 ) by Monte Carlo, it suffices to simulate independent 



replicas of the random variable 3^ := F{Xi)N{X), where hereafter we denote 



N-l 



NiX) =Y[{1 p{Xt, , ^T.+i , r.+i - T,)) ■ 



The exit probability p{x, y, t) does not typically admit a closed form expression and some type of approximation must 
be applied for its evaluation. The short-time asymptotics ( 2.8 ) yields the following natural estimate for p{x, y, t) when 
x,y e D: 



p{x,y,t) := {p{x,y,t)yO) Al, with p{x,y,t) := 



t^ 



s{u — x)s{y — u) 



2 Jia,br ft{y-x) 



du. 



(4.6) 



We also set p{x, y,t) = I ii x ^ D or y ^ D. This approximation satisfies 

P{x,y,t)-p{x,y,t) <ep{x,y,t), (4.7) 



where the error bound ep{x, y, t) is defined as in Remark 3.3 for x,y £ D and by ep{x, y,t)— Oii x^D or y^D. We 
can then approximate N(X) by 

N-l 

N{X) := n (1 - ^(^^^ ' ^^.+1 ' ^^+1 - ^')) ■ (4.8) 

2=0 

Replacing the true exit probability p(a;, y, t) with its approximation p{x, y, t) introduces a bias into the evaluation 
of N{X), which is hard to quantify if the process X is discretized using the uniformly spaced grid T^ = i/N . For this 
reason, we now propose an adaptive algorithm for the determination of the sampling times, which starts by simulating 
the terminal value Xi and then refines the sampling grid, using more discretization points when the estimate of 
the approximation error is "large" . The algorithm is parameterized by a real number 7 > 0, which represents the 
error tolerance and ensures that under suitable conditions on Cp, the global discretization error for approximating the 



quantity of interest (4.1) will be bounded by 7 (see Proposition 4.2 below). The algorithm also requires simulation 



from the marginal distribution /i of Xi and the bridge distribution of X1/2 conditioned io Xt — y (t > 0) . Hereafter, 
we denote the density of this bridge distribution by /f/2(2;, y) and recall the following well-known formula: 

fU-^y) := ^-^^^^^"^^^^ (4.9) 

At the end of this section, we introduce a new method to simulate variates from the density ( |4.9[ ). 

The procedure to generate the skeleton of X is outlined in pseudo-code in Algorithm [T] below. Assume that this 



algorithm terminates in finite time a.s. (see Proposition 4.2 for sufficient conditions for this to hold). The algorithm 
then defines a pair N and T '■— (Tq, . . . jTn), which satisfies the conditions of Lemma |4.1[ Indeed, by construction, 
each Ti takes values in the dyadic grid {i2^™, i — 0, . . . , 2™, m — 0, 1, . . . }, which is a countable set. To check the 
second condition of the lemma, we fix n and a partition tt := {sq, . . . , Sn} of [0, 1], and proceed as follows to write the 
event E := {N = n,To = sq, 21 = Si, . . . , r„ = s„} in terms of {Xs.}"^g: 

• We can and will assume with no loss of generality that tt is a recursive dyadic partition, meaning that {0, 1} C tt 
and, for every t G (0, 1) n tt, there exists fc € N with 2'^t G N, and if we take the smallest such k then also 
i + 2T7 G TT and i — ^ G vr. By construction, if tt does not have this property, the event E has zero probability. 

• We shall assume that n > 2 because if n = 1 then necessarily sq = and si = 1 and, therefore, 

E = {X,^D}U {Xi G D, CpiXo, Xi, 1) < 7} e ^(^0, ^i)- 

• For each £ G {0, . . . ,n — 1}, define tti := {si G tt : 2"~^Si is an even integer}. The number of elements of iTe is 
denoted n^ and the sorted elements of tt^ are denoted s{ < ■ ■ ■ < s^^ . Clearly, ttq = tt and tTu-i 7^ tt since 1/2 G tt 
whenever n > 2; we let £* = maxjZ > : tt/ = tt} and tt* = tt\ tt^.+i. 

• For each i = 1, . . . , n^ — 1, define the event 

{u::ep(xM,X,.^Juj),sl^,-sl)<j{si+,-si)}, if ^ n (sf , sf^^) = 0, 
{oj:eplxAu:),X,.Juj),si^,-sij>j{si^,-si)}, if ^ n (sf, sf+J ^ 0. 

Then, it follows that 

{n n—1 rii — l ^ ( n—1 ni — l 

f]{x,^eD}nf] f]Enul\J{x,,(^D}n fj {x^eDjn fj fj ^' 
1=0 e=f i=l ) [sGtt* sETTt.+i l=e.' + l i=l 

which clearly belongs to cr{Xs^ : i — Q, . . . ,n). 
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Algorithm 1 [X] — GenerateSkeleton(7) 
iVo = 0, A^i = 1, m = 1 
TqI = 0, Ti = 1, Xo = 

Generate an observation Xi from the density /i 
while N„i 7^ iVm_i and {Xt;" G £>, for i = 1, ...,Nm} do 
n = 0, To"+i = 
for i = ^ iVm — 1 do 
AT = 7^™! - T™ 
if ep(XT™,XT™^, AT) > 7AT then 

rpm+l (rpm 1^ T^m \ /rj rpm+1 rprn 

^n+1 — (-'-i "r^i+lJ/^i -'- n+2 •— ^i+l 

Generate an observation Xj,m+i from the bridge density f^x/2^'^ -^t:^i ^ ^t^) 
n ^ n + 2 
else 

-'-71 + 1 ■— ^i+1 

n = n + 1 
end if 
end for 

N,n+i = n 

m = m + 1 
end while 
RETURN X = {(7;'",Xt™)}£B- 

To see that Xj,m+i can be sampled from the bridge density f^rp,^{-,XT'^ — Xt™) in Algorithm 111 we can apply the 



second part of Lemma 4.1 to the couple (fc, 7fc), where Tfe = {Tg, . . . , T^.} contains the first k + 1 sampling times which 
have been added to the grid by the algorithm, in increasing order. 

Algorithm IT] terminates when at least one of the sampling observations Xxi is out of the domain D or the error 
over each subinterval of the sampling times = Tq < • • • < T/v = 1 is small enough in the following sense: 

ep{XT,,XT,^,,T,+i-%)<j{T,+i^T,), i^O,...,N-l. (4.10) 

At first glance, it is not obvious that the algorithm will actually terminate in finite time. The following result gives 
conditions under which this is the case and shows that the global error of the estimate is of order 7. 

Proposition 4.2. The following assertions hold: 

(i) Let X be a Levy process satisfying one of the following two (non mutually exclusive) conditions: 

1. X does not hit points; that is, ¥{t^^' < cx)) = for all x, where t^^^ := inf{s > : Xg = x} or, equivalently, 

3? -y^ ] du — CO, 

where ip{u) = logE[e*"'^i] (see f2^, Theorem 7.12]); 

2. X is a finite variation process. 

Additionally, assume that the upper hound of the approximation error ep{x,y,t) satisfies 

lim- sup ep{x,y,t) = 0, ^a' ,b' e {a,b). (4-11) 

tiO t x,ye(a',b') 

Then, Algorithm[^ terminates infinite time a.s. 
(ii) Assume that K\F{Xi)\ < cx). Let X — {(T^, Ay. )}^q be a skeleton of X on [0, 1] satisfying (4.101 andM{X) be 



given by (4.8). Then, 

|E[i^(Ai)l,>i] - E[F(Ai) N{X)]\ < 7E[|F(Ai)|]. (4.12) 

11 



Remark 4.3. In view of Proposition \4^ ^[F{Xi)l.,-yi] can be approximated by the Monte Carlo estimator 

1 



M 



k=l 



W)^(;t-W) 



where X^ ' are independent copies of the process X and N{X^^) are corresponding values computed with formula 
This estimator has a statistical error which can be estimated in the usual way, and a discretization bias, which 



is bounded from above by 7E[|i^(Xi)|]. In view of (4.13) below, a more precise a posteriori estimate of the bias is 



M 



M 



k=l 



.(k) 



Ak) rr,(k) 



^ Y: l^(^f ^)I 1s(^o E ^p {^Tt^ ' ^rk ' ^+^1 " ^^ 



(K 



with SN:={{XT„,...,XT^)eD''+^} 



Lemma 4.4. Let X be a Levy process such that for all t > Q, the law of Xt has no atom. Then, for all a; € M, 

P[{t e [0, 1] : AXt ^ 0, Xt- = x} = 0] == 1; P[{t e [0, 1] : AXt ^ 0,Xt = x} ^ id] ^ 1. 

Proof. We only prove the first identity, the second one follows by similar arguments (or alternatively by time reversal). 
Let Nl = #{i e [0, 1] : \AXt\ > e, Xt- = x}. Then 



P[{t e [0, 1] : AXt + 0,Xt_ - x} ^ 0] < P\Xl\ < Y,^Wh 
But by the compensation formula (see ^ section 0.5]), 



n=l 



E[iVf] =E 



J\y\>E 



lx,=xv{dy)ds 



v{dy) I V[Xs = x\ds = 0. 

\v\>e Jo 



n 



Proof of Proposition 4-2 Part (i). With the aim of obtaining a contradiction, assume that the statement of 
the proposition is not true, and the algorithm does not terminate. Let {Ti\i>i be the infinite sequence of different 
sampling times produced by the algorithm (in the order in which they were generated, that is, not necessarily ordered 
in time). Let Xi :— Xf be the corresponding sampling observations. Since the sequence {Ti\ is bounded, we can 
find indices {ik}k>i such that Ti^ — )■ T* . Moreover, since every point Ti (for i > 2) is obtained as a midpoint of a 
certain interval, we can find two sequences {T~} and {7;+} such that T:^ t T*, T+ | T*, T* G [T:^ ^T+] for all i and 
and ep{Xj,- ,Xj,+ ,T^ — T^) > "f{T^ — T^) for all i. In addition, since the process X has right and left limits, both 



\miXr, 



X^ and limXrr 



X exist. There are three possibilities. 



li X e (— oo, a) U (6, oo) or X^ e (— oo, a) U (6, oo) then for some i, Xi ^ D, so that the algorithm must have 
stopped in finite time and we have a contradiction. 



li X G (a, 6) and X+ G {a,b) then, using the property (4.111, we can find a contradiction with ep(Xj,- ,Xj,+ ,rj — 

T-)>^{T+~Tr). 

It remains to treat the case when X^ or X+, or both, are at the boundary of D. Then, either X^ — X^ — Xt- 
or AXt' ^ 0. The latter case is ruled out by Lemma [4. 4| and in the case when X cannot hit points, the former case 
is ruled out as well. 

We may therefore assume that X is a finite variation process with nonzero drift // (cf. J9, Theorem 7.12]) and, to 
fix the notation, that X^ = X^ = Xt* = b. We may also assume that T* is irrational, since for every t G Q n [0, 1], 
F[Xt = b] = 0. The fact that T* ^ Q implies that T~ < T* < T^ for every «, and we can also assume that Xrp+ and 
Xj,- belong to D for each i, because otherwise the algorithm would have stopped in finite time. 

Introduce two sequences of stopping times: 

To = inf{i > : Xi > 6} A 1, ct„ = inf{t > r„ : Xj < 6} A 1, t„+i = inf{f > ct„ : Xt > 6} A 1, n> 0. 
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The sequences {t„} and {(T„} do not have an accumulation point except t = 1 and for each n > 0, cr„ > r„ if t„ < 1 
and Tn+i > <Jn if <^n < 1- This holds because for a finite variation process X with drift ^ 7^ 0, {0} is irregular for [0, 00) 
if /i, < and for (—00, 0] if /u > [IHl Theorem 43.20], and X may only creep in the direction opposite to the drift [2^ 
Theorem 7.11]. Then clearly, for every r e [0, 1] such that Xr — b, either there is n > with cr„ — t, which means 
that for some e > 0, Xt ^ D ioi t G {t — e,T), or there is n > with r„ = r, which means that for some e > 0, Xt ^ D 
for t G {t,t + e). In both cases, there is a contradiction with the fact that Xrp+ and Xj,- belong to D for each i. 

Part (ii). Below, we denote p{x,y,t) :— 1 — p{x, y,t), p{x, y,t) — 1 — p{x, y, t), and Sn '■= {(Xtq, ■ ■ ■ t^Tn) G D^^^}. 
Then, since 

N{X)-N(X)^ \{p{Xt^,Xt^^,,T,+^-T,)^ []|(XT.,XT,^,,r,+i-T,), 



we get 



N-l 

|E[F(Xi)l,>i] -E[F(Xi)iV(A')]| <E 



(4.13) 



|F(Xi)|l5„ Y. ep{XT^,XT,^„T,+,-T,) , 

i=0 

which can be bounded by jE [\F{Xi)\]. D 

Simulation of Levy bridges The adaptive method presented in this section requires fast simulation from the bridge 



distribution of Xt/2 conditioned to Xt = y (with t > 0), whose density is given by (4.9). We now propose a simple 
yet efficient method for simulating from the bridge distribution, valid for Levy processes with unimodal density at all 
times. As remarked in Section |3j a sufficient condition for a Levy process to have a unimodal density for alH > is 
that it belongs to the class of self-decomposable processes which includes most of the parametric models used in the 
literature. The algorithm is based on the following simple estimate. 



Proposition 4.5. Let X be a Levy process such that the density ft of Xt is unimodal for all t > 0. Then, 

My) 



ft/2i^, y) < ^-^^YTV^ max{/,/2(x), ft,2{y - x)}. (4.14) 



Proof. For all x and y, 

ft/2{x)ft/2{y ~x) = max{/t/2(x), ft/2{y - x)} mm{ft/2{x),ft/2{y - x)} 

By the assumption of unimodality, the density ft may not have a local minimum, hence, for all a, b, min(/t /2(fl)i ft/2{b)) < 
ft/2 (^) and the result follows. D 

As a consequence of the previous result, random variates with density ftj2{x, y) can be simulated using the classical 
rejection method [T3|, with the proposal density given by f{x) — ^{ft/2{x) + ft/2{y ~ x)), provided that the following 
two requirements are met: 

(a) random variates with density ft{x) can be simulated in bounded time; 

(b) the density ft{x) is known explicitly or can be evaluated in bounded time. 

Assumptions (a) and (b) are satisfied, e.g., for the variance gamma process, normal inverse gaussian process, or for 
stable processes. Simulating a random variable X with density f{x) — ^(ft/2{x) + ft/2{y ~ x)) is straightforward: 
simulate a random variate Z with density ft/2 and an independent Bernoulli random variate U; then, take X = Z ii 
U = and X = y — Z otherwise. 

The expected number of iterations needed until the acceptance for a given value of y is equal to C = <■ ( / ■ 
This number is bounded for Levy processes with Pareto tails such as stable. For processes with lighter tails it may 
be unbounded for large y, but the probability of having a large value of y in an adaptive simulation is very small. 
For example, if we want to simulate Xt/2 and Xt by first simulating Xt and then Xt/2 from the bridge law using 
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formula (4.14), we find that the conditional expectation of tlie number of iterations given Xt equals /(x) ' ^^'^ 



the unconditional expectation is 

-2/,/2(Xt/2) 



E 



5 Numerical illustrations 



ftiXt) 



-2 / ft/2{x/2)dx^4. 



In this section, to simplify the discussion, we assume that the interval D is of the form D = (— oo, b). For the numerical 
implementation of Algorithmlllgiven in Section|4l one needs to be able to perform the following computations efficiently: 



• 



• 



Simulation of the increments of Xt for arbitrary i; 
Evaluation of the density /( of Xf for arbitrary t; 

• Evaluation of the "incomplete convolution" of the Levy density: C{b, y) :— ^^ s{v)s{y — v)dv; 

• Evaluation of the error bound ep{x,y,t), appearing in Algorithmll] 

These computations can be performed relatively easily, for example, for a-stable Levy processes with Levy density 
s{x) = |x|~"~^(c_l3,<o +c+la;>o) and for the variance gamma process with Levy density s(a;) = |x|~^(ce~'^-l^lla;<o + 
cg-^+kll^^p), Por a-stable processes, the increments can be simulated with an explicit algorithm (cf. [Hj), the density 
can be computed using a rapidly convergent series [35] or expressed via special functions (cf. [13]), tabulated for i = 1 
and computed by the scaling property for other values of t. The incomplete convolution is given by 



C{b, y) = c+c^b-'-^^Bil + 2a, 1)F fl + a, 1 + 2a, 2 + 2a, ^1 , (5.1) 

where B is the beta function and F is the hypergeometric function, for which a rapidly converging series is available 
[23] and which can also be tabulated prior to the Monte Carlo computation. For the variance gamma process, the 
density is explicit and the increments are straightforward to simulate [12j . The incomplete convolution is given by 

C(b,y) = — {e-y^+Ei(X(b - y)) - ey^-EUXb)] 

y 

where Ei(x) :— J - — dz, which can also be tabulated, and A := A_ + A+. The error bound Cp for the a-stable or 
the variance gamma process can be obtained along the lines of the general computation of Section [3] or the specific 
computation for the Cauchy process in the Appendix [C] 

For the numerical simulations in this section we shall concentrate on the Cauchy process, which is an a-stable 



process with c+ = c_ := c and a = 1. For this process, formula (5.1) simplifies to 

1 /y\"l c2 (^ b 26^ y 26^ 

n + 
n=l 



C(M) = ^ 1 + 3^;^© =6^|l + ^ + ^ + 6^ + ^l°Kl-f)} 



Note that for small y, the series representation has more stable behavior than the exact formula. The error estimate 
Cp is computed as explained in section [Oof the Appendix. In both examples below, we take c = 1. 

Example 1. In our first example, we evaluate the probability P[supQ<g<]^ AT^ < 1] = P(r > 1), which can be expressed 



in terms of the function (4.1) by taking T = 1, F{Xi) = 1, and the domain (a, 6) = (— oo, 1). Note that in this case, 
the starting value of the process is relatively far from the boundary, and hence the advantage of using the adaptive 
algorithm is less important. The process will typically cross the boundary by a large jump with a large overshoot, 
which makes the exit easy to detect, even with a uniform discretization. 

We study the performance of our adaptive algorithm for various values of 7, and compare it to the standard uniform 
discretization. When interpreting the results of simulations, one needs to distinguish between the actual error (that 
is, the difference between the computed value and the true value), and the theoretical value of the bias (computed as 
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Time 



Figure 1: Illustration for Example 1. Left: values returned by the uniform discretization algorithm and the adaptive 
algorithm, as fonction of the computational time for 10^ paths, measured in seconds. Different points on the graph 
correspond to different numbers of discretization dates for the uniform discretization (ranging from 32 to 8192) and 
different values of the tolerance parameter 7 for the adaptive algorithm (ranging from 7 to 7 x 10"'*). The curve for 
the uniform discretization is smooth because all the points have been generated using the same trajectories, while for 
the adaptive discretization different paths have been used. Right: comparison of the theoretical bias of the adaptive 
algorithm with the actual discretization bias of the uniform discretization. 



explained in Remark 4.3 above), which does not require the knowledge of the true value. As an estimate of the true 



value, we use the value computed in an independent simulation by uniform discretization with 16384 points and 10^ 
trajectories, which is approximately equal to 0.38935 with a standard deviation of lO""*. The difference between the 
values for 8192 and 16384 points (on the same trajectories) is smaller than 10~^, hence one can presume that, for all 
practical purposes, convergence up to this precision has been achieved. 

Figure [1] shows the dependence of the values computed by the two algorithms on the computational time required 
for 10^ MC trajectories, for different numbers of discretization points (for the uniform discretization) and different 
values of the tolerance parameter 7 (for the adaptive algorithm). While the uniform discretization algorithm exibits 
a clear bias which decreases as the number of discretization dates increases, the adaptive algorithm removes the bias 
completely; all values returned by this algorithm are within confidence bounds of the true value. 

The theoretical bias, computed as explained in Remark |4.3[ is greater than the actual error, because the error 
estimates of Appendix [C] are upper bounds, and because it does not take into account the possible cancellation of errors 
on different intervals. Figure 111 right graph, compares the theoretical estimate of the bias of the adaptive algorithm 
with the actual bias of the uniform discretization. One can see that for small computational times, the theoretical 
bias for the adaptive algorithm is greater than the error of the uniform discretization, however, the theoretical bias 
converges to zero much faster, and for relatively large computational times is actually smaller than the error of the 
uniform discretization. The empirical convergence rate (estimated from the slope of the straight lines) is T~°-®* for 
the uniform discretization and T^^'^ for the theoretical bias of the adaptive algorithm. 



Example 2. In our se cond example, we evaluate the probability P[supo<5<i ^s < 10 ^] , which again can be expressed 
in terms of the function (4.1 ) by taking T = 1, F{Xi) = 1, and the domain (a, b) — (—00, 10~^). In contrast to Example 



1, here we consider a situation where the starting point is close to the boundary. In this case, as we shall see below, the 
advantage of the adaptive algorithm is more striking, since the process can cross the boundary and come back while 
it is still close to the starting point and, hence, a very fine discretization will be necessary to detect this event with 
uniformly spaced observations. As a result, for the uniform discretization we do not observe convergence to a sufficient 
precision even with 16384 points, and therefore the true value cannot be estimated as in the previous example. Instead, 
we shall use as the true value the value produced by the adaptive algorithm with 10^ Monte Carlo paths and equal to 
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Figure 2: Illustration for Example 2. Left: values returned by the uniform discretization algorithm and the adaptive 
algorithm, as function of the computational time for 10^ paths, measured in seconds. Different points on the graph 
correspond to different numbers of discretization dates for the uniform discretization (ranging from 256 to 16384) and 
different values of the tolerance parameter 7 for the adaptive algorithm (ranging from 9 to 9x 10^'^). Right: comparison 
of the theoretical bias of the adaptive algorithm with the actual discretization bias of the uniform discretization. 

0.0360, with standard deviation of 6 x 10~^ and theoretical bias of 3 x 10~^. 

Similarly to the previous example, Figure [2] shows the dependence of the values computed by the two algorithms 
on the computational time required for 10^ MC trajectories. Here, the adaptive algorithm exibits the same kind of 
behavior as in the Example 1 above: all the points generated by the algorithm are within the confidence bounds of 
the true value. However, for the uniform discretization, the convergence is much slower than before and only the last 
value obtained with 16384 discretization points falls within the confidence bounds. Figure [21 right graph, compares 
the theoretical estimate of the bias of the adaptive algorithm with the actual bias of the uniform discretization. Once 
again, the behavior of the adaptive algorithm is roughly the same as in the previous example, showing that the method 
is robust with respect to the parameters on the problem. On the other hand, as expected, the uniform discretization 
presents a significant bias in this case (the convergence rates are similar to those obtained in the previous example, 
but the constant for the uniform discretization is much bigger). 

Acknowledgments: The authors are grateful to an anonymous referee and both the associate editor and the editor- 
in-chief for their constructive and insightful comments that greatly helped to improve the paper. 



A Proofs of Section [2] 
A.l Proof of Theorem [272] 



Throughout the proof, we shall use the notation 

Yt := sup Ys 

0<s<t 



and Yf 



:= inf n, 

0<s<t 



(A.l) 



for a given cadlag process (lt)t>o- Without loss of generality (by considering separately the positive and the negative 
part), we can and will assume that ip is nonnegative. Additionally, assume that a £ (—00, 0) and b e (0, cx)). The cases 
a — —00 and 6 = 00 will be evident from the proof below. We also let Htpljoo '■= ess snp^ip{x), ||iy9||Lip be the Lipschitz 



S,y + d), J] := <5o/2. 



sup,. 



iXf 



norm of ip, Is{y) := {y 

and flg :— sup^. s^{x), which are finite in light of (2.1). In what follows, J^ '■— '^ i-^s ■ s <t)\/ J\f where Af denotes the 
null sets of T. To lighten the notation below, whenever the ess sup of a function g, defined £-a.e. in some region, is 
considered, we shall simply write sup„5(u) instead of ess sup^g(u). 



b A |a|, B := {t <t} = {Xt > 6 or A^ < a}, f/f 
In what follows, J^ '■— '^ i-^s 



16 



The idea is to condition on the number of jumps of the compound Poisson component Z^ . To this end, let us denote 

Ak{t)^'&{Lp(Xr)\{r<t,Xt<:ill,lyy').Nl=k)) , for fc = 0, 1, 2, and A-i{t) =^{ip{Xr)\{r<t,Xt<^Ii,{y),Nl>?,-}) ■, 

so that clearly 



E((^(X,)l{,<t,x,e/,fe)}) = ^0 W + • • • + M{t). 
Note that each of the terms on the right-hand side of the previous equation can be expressed as 

Akit)^ \ P^{u)du, (fc = 0,...,3), 

Jy-& 



{A.2) 



(A.3) 



for some nonnegative functions Pi[u). Indeed, for fc = 0, 1, 2, by the standard definition of conditional expectation, 






(A.4) 



ly-S Jy-S 

The case fc = 3 is treated in the same way. Let us analyze each of the four terms in the right-hand side of (A.2[). 



(1) No big jump. Note that, on the event N^ — 0, Xg — X^ for all s < i and, thus, {t < t} = {r'^ < t}, where 
T^ := inf{u >0:X^(^{a, b)}. Therefore, 

Ao{t) = E (^(X^Ol{r^<*,X|e/,fe),iVf=0}) - E {ip{X'r^)l{r^<t,X^^I,(y)})P (^ = 0) , 

where in the last equality we used the independence of X^ and N'^. Next, conditioning on T^e, it follows that 
By Markov's property, 

where F{z, s) = P (z + ATf e hiy))- Note that if t^ = t, then F {X^,,t - r'^) = since X^, e (a, 6)^= and Is{y) C (a, h). 
On the other hand, on the event r^ < i. 



,y+5 



rV+S 



F{X^,,t-T')^ f^_^,{u-X^,)du< sup sup fl{u~x)du, 

Jy-S Jy-S 0<s<t xeia.b)" 



since again X^^ G {a,by. Putting the two previous cases together and recalling (A.3), we have 



Aoit)= I P°{u)du< I (e-^=*||(p||oo sup sup fl{u-x)]du 



y-S 



y-6 



0<s<txe(a,b)'= 



y+S 



Pt{u)du, (A.5) 



y-d 



implying that Pi{u) < P^{u), for £-a.e. u E {a + SQ,b ~ Sq). Furthermore, using (2.4 



sup Pt{u)< sup P°(u)< ||</J|Uc3((5o,e)t^ (t < io)- 

a+So<u<b—So a+Sa<u<b—5a 



(2) One big jump. Let r^ and Yi be the time and size of the i jump of Z^ . Clearly, on the event {N^ = 1} 

^{XT)'i-{T<t,Xteh{y),N^ = l} = V{^T)'^{T<Ti,XI+YieIs{y).N{ = l} + ^i^r + Yl)l{Ti<T<t,X^+Yiehiy),N^ = l} 

< \\^\\,x,'i-{X^+Yiel6{y),N^ = l} \^{Xf>boi- Xl<a} + '^{X^+Yi>b or Xt+Yi<a} ^ 



17 



It follows that 



< Ai{t) < \\(p\\^E {'i-{Uf>c,XI+YieIs{y),N- = l}) + ll</'ll<x>E \^l^Yi>b-X^ or Y^<a-Xl}'^{X^+YieIs{y).N- = l} 

= er^^'^WfW^X.tE {l{u,>c,XI+Yieh{y)}\+ C'^^'*y\\^^etK [l^Y^>b-X^ or Y,<a-X'^}'^{XI+Yieh{y)} 



^i,i(*) 



Ai,2{t) 



where in the last equality we use the joint independence of A^^, Yi, and X^ . Conditioning on a{Xg : s > 0) and 
applying Fubini, 

^i.iW^e-^-l'^llooiE l{t/f>c} / s,{v)dv]^ e-^^^%^\\^tK{l{ui>c}Se{u-XI))du. (A.6) 

V Jy-S-Xf J Jy-S -^ ^ ' 

Pl-\u) 



Using (2.1) and Lemma 2.1 sup„P/'^(u) < e~^^^t\\ip\\ryoaeV {U^ > c) < e''^^^a^\\ip\\aoC2{c,e)t^ , where e > is chosen 
small enough. Similarly, conditioning on cr(X| : s > 0), making the substitution u = X^ + w, and applying Fubini, 



(A.7) 



= / e^^''^\W\\ootV.[l{^^a+x;-X'; oru>b+Xi-X'r}Se{u- Xl)\ du. 
Jy-s V 1 .. i. 



Pl'\u) 



Using again Lemma |2.1 



sup Pt'\u) < e-^-*||(p||ootoeP [X! -Xl>5o or X^ - X^ > 6o) < e-^''Moota,P Xf - X^ > So 



uG(a+i5o,b— i5o) 



< e-^^'\\ip\\^ta,¥ ( suplXfl > S„/2 ] < e-^^'M^a,C2iSo/2,e)t^ 



(A.8 



Therefore, recalling from (A. 3), the nonnegative function Pl{u) is such that, for £-a.e. u E {a + So,b — do), < Pl{u) < 
ELiPl''i^)<MooaetHC2{c,e) + C2iv,e)). 

(3) Two big jumps. As before, let r^ and Yi be the time and size of the i*'* jump of Z^. Clearly, 

yy{Xr)l{r<t,Xteh{y),N^=2} = 'PiXr)i{r<Ti,X} +Yi+Y2el5iy),Nf=2} 

+ Vi^r + yi)'^{Ti<T<T2,Xf+Yi+Y2eh(y),Ni=2} 
+ ^{Xl +Yi+ Y2)l{r^<r<t^XI+Yi+Y2(iIs(y).NI=2} 
< \W\\oo'^{3s<Ti:X'^(^{a.b);X'i+Yi+Y2el5{y)\NI=2} 

+ Vi^l + ^l)l{3se[Ti,r2):Xf+Fi^(a,6);Xf+Yi+y2e/a(a);Aff=2}' 
+ \\'P\\-xi'^{5se[T2A-Xl+Yi+Y2i(a,by,X^+Yi+Y2eh(y)-,N^=2}- 

Then, using the independence of N"^ , the l^i's, and X'^ in the first and last terms, we have the inequality: 
A2{t) < e-^^'{ty2)\lM^E{l^u;^>,^xi+n+Y2eis{y)}) 

+ E yp{X^ + yL)l{Xf +Yi>b or XI+Yi<a;X^+Yi+Y2eIs{y);N^=2}j (A-9) 

+ e " (i /2)Ag||(y9||ooE (^l{X|+yi+y2>6or Xf+yi+i'2<a;-^f+i'i+y2e/«(y)}J' 
=: A2a(i)+A2.2(t)+A2^3W■ 
As before, conditioning on cr{X^ : s > 0), changing variable from w to u ~ X^ + w + w, and applying Fubini, 



(A.IO) 



A2,l{t) = e - 2 ||<^||ooi I^ ( / / '^{Ui>c}'^{y-S<Xt+w+v<y+S}Se{v)Se{w)dvdw 

y-\-5 noo ry-\-S 

e-^-*2-i||^||ooi' / s,{v)E{l{u.>,}S,{u-XI-v))dvdu=: / Pt'\u)du, 

y—S '^— oo J y—6 
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and, hence, 



sup Pt'iu) < e-^''2-'Moot'XeaePiU^ >c)< e-^''2-'M^X,aeCi{c,e)t 



Similarly, 



^2,3(0 = e "2 yWoot ^[ ^Xl+v+w>b oi-X',+v+w<a}^{v-S<Xi+iu+v<y+S}Se{v)s^{w)dvdw 

= / e-^^*2-^\\ip\\oct^ / E (l{x}-xf+u>b or xi-xi+u<a}S6{u - Xf - w)j s,{v)dvdu 



y-S 
y+S 



Pt'^iu)du, 



(A.ll) 



and, thus, as in (A 



sup Pt^iu) < e-^''2-'M^t'X,a,V {XI - Xt > So or X^ - X^ > So) < e-^^*2-'MooXeaeCiiSo/2,e)t\ 

u£[a+5a,b—Sa] 

Finally, we provide an upper bound for A2,2it)- First, we use the bound (p{X^ + Yi) < (p{Yi) + \\ip\\upUf and again 
the independence of N'^, the Yi's, and X^ to get 

^2.2 (t) < e-^^*(tV2)A2E ({^(Fi) + MlupUn l{Xf+yi>6 or X^,+Y,<a-Xl+Y,+Y,els{y)}) ■ 

Next, by conditioning on a{X^ : s > 0) V (j{Yi), we may writer] 

/ .y+s-xl^Yi \ 

A2At) < e-^»*(iV2)AeE {^(Yi) + M\upUn l{xf+n>fc or x',+Y,<a} / s,{w)dw 

( r .y+S- XI -V \ 

= e-^'*(iV2)E| / {v{v) + Mu^UI}s,{v) I s,{w)dwdv\ (A.12) 



(a-X-A-X-Y 



y-S-X^-v 



y-S 



g-A,t2-lt2g 



(a-Xl,b-XlY 



{^{v) + lltysllLipC/t"} s^{v)se{u - X^ - v)dv du. 



P?'^H 



In order to find a lower bound for ^2(^)1 note that 

y^iXr)l{T<t,Xt€ls{y),Nl=2} > ^{X^ +Yi)l[ri<T<r2,X^+Yi+Y2ehiv):N^='2} 

> fiX^ + Yi)l^Yi+XI>b or XI+Yi<a}^{X^<b,X't>a}^{XI+Yi+Y2ehiy),Nl=2}- 

Using the previous inequality and the lower bound fiX^ + Yi) > (p{Yi) — \\ip\\hipUf together with the independence of 
N", the Yi's, and X"", it follows that 



^2(i) > e 



-A.t(Aei)^ 



y+S 



E [{f(Yi) - \\(p\\LipUt}l[Yte{a~XfM-X^^)'=,Xf<b,XI>a,XI+Yt+Y2eh{y)}) 

{ip{v) - ll^llLipC/f } s,{v)s,{u - XI - v)dv ] du. 



y-S \ J{a-X^,b-Xl) 



P?H 



2,2/ 



As it will be proved in Lemma A.l below, P^ ' (u) and P^ (u) arc such that 



lim sup 



lim sup 



1 d2,2/ X 1 

-^ P, (u) 



ip{v)Si.{v)s^{u — v)dv 



{a,bY 



t^o 



ue{a+5o,b—So) 



1 1 /■ 

T^P^iu)-^ ^p{v)s^{v)s^{u-v)dv 

t 2 J(„ f,)e 



= 0. 



(A.13) 

(A.14) 



^Here and below we use the convention {x, y) = and {x, y)'^ = {—00, 00) for x > y. 
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Using (A.9), (A. 14 1 and the corresponding bounds for P^ ' (u) and P^ ' (u), it follows that the nonnegative function 



Pt{u) defined in (A.3I is such that 



lim sup 

*^Oue(a+<5o,b-<5o) 



1 1 /■ 

T^Ptiu) - ^ / fiv)seiv)s,{u - v)dv 



= 0. 



(A.15) 



(4) Three or more big jumps. As before, we have the following bound 

= ||</j||ooP(iVf =n) [ E (sr(w - X[)) du, 

Jy-S 



Using the previous inequality and (A.3), we have 

^fy+s 
y-s 



ry+S ry+S r °° ^n 



■y+s 
du =: / Pf(u)du. 

y-S 



Since Prlloo < Ar'as 



i" 



supp3(«) < e-^^*a,||^|U E -K'-' < C{e)t^ 

n— J 



(A.16) 



for some constant C(e) < oo, and we conclude that < Pfiu) < C(£)i for £-a.e. u 

Putting the foui 
Pt{u) such that 



Putting the four previous steps together, we conclude that E {ip{Xr)'i.{r<t,Xte(y-s,y+S)}) ~ /^-a Pt{u)du, for a function 



lim sup 



t->o 



ue{a+SoM—So) 



^Pdu)-l 



{a,by 



(^(w)Se(w)Se(M — v)dv 



= 0. 



Finally, it is easy to see that for any m e (o + i5o, ^ — 6o) and a < < 6, there exists an Eq > small enough such that 
/(o b)" v{''^)^e{'>^)seiu — v)dv = /, j^x^ (p{v)s{v)s{u — v)dv, for all < e < Eq. This concludes the proof. D 

Lemma A.l. Verification of {A. 13) and (A.IJ;.). 



Proof. Let < e < 1 and M^ :— Xf — /i^i be the martingale component of X"^ . We shall analyze the expressions 



appearing inside the absolute values on the right hand side of equations (A. 13) and (A. 14). Define the random intervals 



/ := (a — Xf , h — Xf), / := (a — Xf ^ b — 2Lt), and the corresponding limiting interval J = (a, b), under the convention 
(x, y) ~ 9 ii y < X. Denote 



DUu)^e(J y{v) + \MupU!}se{v)s,{u~Xt~v)dv\ ~ J ^{v)s,{v)s,{u~v)dv, 
D'^{u) ^E(l^x-<b,x-[>a} IW{v) - \\Lp\\upU!} Se{v)se{u ~ XI - v)dvj - / (p{v)se{v)s^{u-v)dv. 
Let us first analyze Dj. Clearly, 

Dl{u)^\\ip\\upK(ui s,{v)se{u-XI -v)dv] +e( / (p{v)se{v)seiu - X^ - v)dv 



+ E( I ip{v)si,{v)[se{u - Xf - v) - Ss{u - v)]dv ] , 



and, therefore, using that P \ J'' d (a, a ~ X_f) U (6 — X^ , 6), under the convention (— oo, — cxi) = (oo, oo) = 
\Dl{u)\ < a,A,||^||LipEC/,- + alMUEiX^ - XI) + X,y\U\s'JUE\XI\ 

< (aeAe||(^||Lip + 2M^al){Esup \MI\ + \^i,\t) + ||(^||ooAe||4||oo (E|Aff I + |/ie|i) 



s<t 
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Using the trivial inequalities (Esup^<t |A//J^|)2 < Esvip^^^{M|f and (E|AfJ^|)2 < E(A/J)^ together with Doob's in- 
equality, we then get the bound 



(A.17) 



\Dl{u)\ < [2a,\,y\\up + MMocal + \\(p\\ocXe\\s',\\oc]<Jet^^^ 
+ [a^X^\\ip\\up + 2||</j||ooae + ||¥'||ooAe|ls^||oo] IfJ-elt, 

where ct^ := <j^ + J c^{x)x'^i^{dx). For Df{u), note that 

Dt{u) = -Ef l{X=>6orXf<a} j y{v)Se{v)s,{u - Xf - v)dv\ + \\f\\-LipM'h-{x- <b,Xl>a}U! j Je{v)s^{u - Xf - v)dv 

+ E / {p{v)s^{v)s^{u ~- Xi ~ v)dv\ — I ip{v)s^{v)ss{u — v)dv. 



Defining c = \a\ Ab and following the same steps as above, it is easy to verify that |_D^(u)| admits the following upper 
bound: 

\dUu)\ < |l^|UAe«.P((7f > c) + a,A,|l(p|lLipE[/f + 2a'j^\\^EU! + A,|l^|Ull4llooE|Xf | 

< \\ip\\oo>^eaeCi{c,e)t+ [2aeAe||(y3||Lip +4||(p||ooa^ + ||(^||ooAe||4l|oo]o-£t^^^ 



+ [aeA^ll^llLip + 2\\ip\\aaal + ||<y5||ooAe||s^||oo] |/^e|t. 



where we had used the tail probability bound in (2.4). 



(A.18) 

D 



A. 2 Proof of Proposition 2.4 



We use the notation introduced at the beginning of Section A.l above and, as before, we assume without loss of 
generality that (p is nonnegative. As it was done in (A.2|, by partitioning the space into the different values that TVf 



can take on, we can decompose E{ip(XT)liXt€h(y)}) into three terms: no big jumps, one big jump, and two or more 
big jumps. These terms can in turn be expressed as integrals of the form (A. 3 1 using a procedure similar to (A.4). The 
term with no big jumps is such that 

y+S py+S 

P?{u)du := E((p(X,)l{x,e/,fe),JV-o}) < lbllooP(Xf e /,(y),7Vf = 0) = / e-^''M^Jl{u)du, 
y-S Jy-S 



■- p->et|! 



|]i^||oo/f (^)- Using (2.4ii), we can further upper 



which yields an upper bound for P°(m) of the form Pt{u) 

bound Ptiu) by ||v'||ooC2(c, e)t'^ uniformly in (a — i5o, ^ + <5o)^- The term with two or more big jumps can be bounded 

similarly to the term with three or more big jumps in the previous section. Concretely, this term satisfies 



py+S 

/ P^{u)du:=E{Lp{Xr)l 

Jy-S 



{Xteis(y),Ni>2}, 



< 



y+S 



y-S 



n=2 



-X.t 



t" 



M^E{sr{u~Xt)) 



y+S _ 

du =: I Pf. {u)du, 

y-S 



and, using that ||s*"||oo < A" ^a^, we can further upper bound Pt{u) by C{e)t^ for a constant C{e) < oo. The term 
with exactly one jump is decomposed as follows: 



y+S 



y-s 



Pl{u)du := E {ip{Xr)l{Xtei.',(y),N-=i}) = E {ip{Xl)l{Xteh(y);r<T^;N-=i}) + E {tp{Xr)l{Xteh(y),T>Tr,N-=i}) 



where ri is the time of the first big jump. Out of these two terms, the first one satisfies 

E(^(X^)l{x,e/,(,);r<r,;iv-i}) < llv'llooPps G [0, i] : Xf i D;Xt + Yi e Is:N! = 1] 

= e-^-*A,i||v?||ooP[3s e [0,t] : X^ i D;X^ + Yi e h] 

< f e-^^Hy\\^E[lui>cSeiu-Xf)]du, 
Jy-s 
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where the integrand P^.' (u) := e ^^*t\\ip\\^E[lifE^cSe{'>J- ~ ^t)] ^^ uniformly bounded by \\ip\\ac,a£Ci{c, e)t'^ . As for 
the second term E {ip{Xr)l{Xteisiy}.r>Tv,Ni=i}) = E {ip{X^ + Yi)l{x-+Yieh{y);T>Ti:N-=i}), it can be bounded from 
above by 

E{ip{Yi)l[X^+Yieh{y);Nl=l}) + \\'P\\hipK{Utl{xi+Ytehiv);Nf=l}) 
l-y+S 

{e'^^HE [(p{u - X't)s,{u - Xf )] + e-^-4||(^||LipE [U^s,{u - Xf )]} du 

ly-S 

< f {t^{u)s,{u) + t{\\ip\\upa, + |l^|U||4l|oo)E[|Xf |] + tMupa,E[U!]} du. 

Jy-S 

Similarly, this can be bounded from below by 

^{'P{Yl)l{x^ +Yi£ls(y);X^ <b:XI>a;N^ = l}) " ll</'l|LipE(C/f l{xf +yie/5(i/);Aff = 1}) 

{ e-^'HEiipiu - X!)s,{u - X[)l^x^<b,xi>a}) - e-^^HMupE{U!s,{u - Xf )) [ du 

ly-S 
py+S 

> / {tLp{u)se{u) - WipWocasKt'^ -^(||<^||Lipa£ + ||(^||oo||sel|oo)E[|^f|] ~ t\\ip\\upaeE[Ut] - t\\ifi\\oca£Ci{c,e)t} du. 

Jy-S 



To conclude, we estimate E[|Xf |] and E[C/f ] as in the proof of Lemma A.l above. D 



B Proofs of Section [3] 



In this part, we provide the building blocks to develop an upper bound for the remainder TZt{u) appearing in (2.5) 



B.l Proof of Lemma 13.11 

Let us first assume that /i > so that Xf := Mt + fit is a submartingale. By Doob's inequality, for all c > 0, 

supX, > 7/^ = P fsupe=^= > e'^"^ < ^}^lJl ^ e'i'{<^)-cn (B.l) 

s<t J \s<t J 6'='' 

with ip{c) = IJLC+ ^^^ + /uKgle*^^ — 1 — cz)v{dz). Minimizing the right-hand side over all c > 0, we get, as in [38] (see 



p. 87 therein) 



inf e*'^(")-"" = exp ( -t / T[z)dz ) = exp ( -i / T{z)dz ) , (B.2) 



'=>o V H'{Q) I \ Jt^ 



where we are taking t < rj/fj, and r : [0, oo) — > M is the inverse function of 

^'{x)= fi + a'^x+ f zie^"" -l)v{dz). (B.3) 



J\z\<e 

As in [25], note that, for a; > 0, 



0< / zie'"" - l)v{dz) < f |z|(el^l"-l)iy(dz)< / |z| V 

J\z\<e J\z\<e J\z\<e ,_, 



{\z\xf 



l\z\<e J\z\<e J\z\<e —^ k\ 



i'{dz) 



< / \z\^v{dz) Y^ ^-^ = / \z\Mdz)- (e- - 1) . 

J\z\<e r-^ *! J\z\<e e 
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From the previous inequality, for x > 0, 



1 



< i}'{x) <fi + cr^x+ / \zfi^{dz)- (e''^ - 1) < M + 

J\z\<s £ £ 

where we used the fact that a"^ = a'^ + J,^,^ \z\'^i'{dz). This imphes that 



e^^ - 1 2 



t(z) > - log 1 + 



z — n 



and therefore, substituting this into (B.l ) and (B.2) and using that wln(w) < (1 + i;) ln(l + v) and e ""log"" < g*^ for 
all w > 0, we have 



supXs > T] 

s<t 



£ .-'0 



tcrf 



log(l + s)ds 
e(?7 - /it) 



ta'i 



log 1 



eiT]-^t)\ s{r]-nt) 



tai 



,exp,-l^,oj£i^)U,?.^-„p 






log 



ecrf 



The above inequality proves the statement (2-i) for the case fj, = 0. Next, it is easy to check that the function 
u — )• {u/s)log{eu/ea'^) is strictly convex in (0,oo) and reaches its global minimum value of — crf/e^ at w = a'^/e. 
Hence, whenever rj — ^t > 0, 



sup Xs > rj 






which proves the statement (1) for /i > 0. Also, if /i > 0, i < f]/^, and rj < a'^/e, we have that 



exp 



V - i^t ^ ( £{V - l^i) 



< sup exp < log 



u . 



e \ ea, 



eu 



-?loK ^ 



ST] 



0<U<7] 

which proves the statement (2-ii). Finally, we consider the case fj, < 0. In that case, obviously, Mt + iJ.t < Mt and 

V 

sup (M^+Ais) > 77) <P (supA/^ > ?y) <i^ ( — ) <i^e>, 

where in the second inequality we used the case (2-i) with /i = that was proved above. The previous inequality proves 
the bounds (2-i) and (1) for /i < 0. D. 



B.2 Proof of Theorem [3:21 



To prove the estimate (3.11 1 for the remainder TZt{y), we analyze each of the four terms in (A. 2 1 contributing to it 



(No big jump) The first component of the error is due to P° which, as seen in (A.5), can be bounded by 



e(°)(0,y,t) :=e ^'^\\lp\\oo sup sup f^{y-x) = e ^=*||(/3||oo sup sup f^{z) 



0<u<t xe{a,b)'= 



0<u<t z^(y-b,y-aY 



Next, recalling the notation Aj, = {b — y)A{y — a) > and employing our hypothesis that Xf has unimodal distribution, 
we can further apply the bound (3.10) to get 



e(°)(0,y,t)<e-^^*^'"'^"'- 



A 



sup 



y 0<u<t 



K\ > ^ 



< '5:^^C(A./4,.),*, 



fort<to(e,Ay/2)Ati(e,Ay/2). 
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(One big jump) There are two sub-components to the error in this case. The first is due to P^'^ in (A. 6). This term 
can be bounded by 

e(i'i)(0, y, t) := |l^|Ue-^=*tE (l{y.>,}S,(y - Xf )) < \M^e-^^Ha,nU! > c) < 2\M^e-^^'a,C[c/2, e)f+^, 

for t < to{e,c/2). The other sub-component is due to P^'^ in (A. 7), which can be bounded, for t < to{e,Ay/2), as 
foUows: 

< M\ooe-^^Ha,¥(snp\Xl\ > ^) < 2M^e-^^' a,C{ Ay /4,e)f+^ . 

\u<t ^ I 



(Three or more big jumps) This component can be bounded as in (A.16): 



e(3)(0,y,t) := ||^|Ue-^^*5^-E(C(«-^n) < MUo.,\-^ (l - e-^^*[l + A,i + (A,t)V2]) • (B.4) 



n— 3 



(Two big jumps) There are three sub-components to the error in this case. From (A.IO) 



e(2'i)(0,2/,i) := MUe^^^'- \ sa«)E {l{c/=>c}Ss(y - ^f - t')} d" < llvlUe~^^*a,A,C(c/2,e)i2+K, 



for t < io(e, c). Similarly, from (A. 11 1 



e(2^3)(o,y,i):^||^||^e 



-A,t' 



,(«)e{] 



-A^t' 



{Xf-Xf+i/>fcorX; 



A, 



%~Xl+y<a}Se{y ~ ^l ~ V)^ 



dv 



< l|(^|Ue-^^^^%AeP(sup|XS| > ^ ) < M\ooe-^^'a,X,C{Ay/4,e)t'+^ 






(B.5) 



(B.6) 



for t < tf){e,Ay/2). Next, we consider the error due to the limits ( A.13||AT4 ). These were bounded in Lemma A.l 
Hence, by taking the maximum of (A. 17) and (A. 18), after some simplification, we get the following expression for the 
error term e(^'^^(0, j/,i): 



Finally, we also need to take into account the error due to approximating e ^^*y // h)<:'p{''j)se{v)se{y — v)dv by 

V I(a bY '/'(^)*£(")se(y ~ v)dv, which is of order \\ip\\ooXla,.t^ /2. Putting all the previous bounds together, we obtain 
the overall bound (3.11). D 



C Finding the estimate ef{0,y,t) for the Cauchy process 

In this paragraph our aim is to find an explicit bound for the Cauchy process with Levy density I'ix) = t^j- (and no 
drift), which is used in the numerical illustrations. For simplicity, we shall only consider the one-sided case (a = — oo). 
Setting Ci;{x) = 1|2,|>£, we get /^^ = for all e, and the law of the process is symmetric, which means that tQ(e,ri) — 



ti{e, rj) — ~\-co for all e > and 77 > 0. Moreover, a^ = 2ce and Lemma 3.1 implies that P[sup^<f -'^t > ??] < i^ C'(r/, e) 
[supj,<j \Xt\ > rj] < 2t^C{rj, e) with C{ri, e) — I ^^ j " . The results of the above section can now be improved to 



and 



4e 



-A,t 



1 + 



e(i'2)(o,y,i)<2||^||ooe-^^*a,C(£,(fe-y)/2)ti+^, e^'^'\Q,y,t) < ^^e-''^'aAeC{b,e)t 



2+ 
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To estimate e'^'^^ more precisely, let £o < '' ^ ^ A (5 — e). Then, 



IA'(y)l < 2||(^IUaeAeP(C/f > eo) + 2||(^||LipE 



HocE 



U! / s,iv - Xt)s,{y ~v + Xt- Xl)dv 



lc/,-<eo / {se{v ' Xl)s,{y ^ V + X^ - X^) - s,iv)s,{y - v))dv 



< 2a,Ae(||(^|looP {U! > So) + |l<^|lLipE[C/t1) - ||(/^||ooIE [C/f ] / s'^{v ~ eo)s,{y - « + 2eo)dv 

b-y. 



+ 2M^E[UnJ s,{v)s'^{y - V + ^)dv 

< 2a,X,{M^¥ {Ut > eo) + y\\upE[U^]) + 2||(p||ooE [U^] s, {b - Eq) s, {b - y - 2eo) ■ 

A similar argument shows that 

\Dfiy)\ < s,{b)X,{2Mupmt] + II^IUP[^t^ > b]) + MUE [C/f] s, (6) s, {b - y) , 

which means that the bound for jZJj (y)| always dominates. Using the former bound, we finally find the following upper 
bound for e^^^(0, y, i): 

2M^e-^''a,X,C{eQ, e)t2+? + 2e-^'Hia,{s, (6 - Eq) Se {b - y - 2£o) llt^lloo + a£A,|l(^|lLip}. 

To speciahze the estimate e'^\ we upper bound A"P(Xt >b,Xt^ Is{y)\N^ — n) by 

^ k n \ n / k n 

A^p(x,- + max Y.Y, > 6, xf + ^ F, e is{t)) < A?^p(x[ + ^y, > &,xf + ^r, e Isit) 



0<k<n ■ 

^ ^ ~ i=l i=l 

The cases fc = and k = n are treated separately: 



fe=0 



1=1 



A^P X^ > b,X^ + J2^^^ ^■'W < -^P supX;; > 6 sups*"(a;) < SC{b,e)t^ sups*"(a;)- 

/ n n \ 

A;T ix|+J2Y^> b, Xf + ^Y,eh{t)\ < SP {XI -XI + y + S>b) sup s*"(a;) < 25C{{b - y)/2, e)i^ sup s*"(a;). 



For < k < n, 



i=l 



^f + ^K, >6,Xf + ^K, G/5(t) =E 



j=i 



dw / dvsl''{v)s;^''-''\u-v-X't) 



<<5sups*"(a;)P(L/f >eo) + <5sf'(^-eo) / dv sl^''-^\y - v + 2eQ + 5) 

X Jb 

where s^ is any function which is increasing on (— oo,0), decreasing on (0,cxd) and satisfies Si;{x) > s^{x) for all x. For 
the Cauchy process one can take 



Se{x) 



2c , ^k, ^ 1 /27rc 

: ^ SO that K (x^ = - I 



ek 



n \ e J [eky + a;^ ' J^ 
Assembling all the estimates together, we finally get 



sf{v)dv 



( 



1 /27rc 



e(3)(0,y,t) < '^^a,\lt\C{b,e)ti +2C{{b - y)/2,e)t 



'^+2C(eo, £)*"?) + 



TT V e 



leTTC^i^lJVS 



arctan ■ 



ek 



-\^t 



3eib-eo)Hb-y-2eo) 



The above estimates satisfy condition (2.6) for e < -^ A b. In the numerical examples discussed in the paper we have 



taken e — —^ A | and Eq 



b-y A b 
4 " 2- 
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